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Interactions between rhenium impurities in silicon are investigated by means of the density func- 
tional theory (DFT) and the DFT+U scheme. All couplings between impurities are ferromagnetic 
except the Re- Re dimers which in the DFT method are nonmagnetic, due to formation of the chemi- 
cal bond supported by substantial relaxation of the geometry. The critical temperature is calculated 
by means of classical Monte Carlo (MC) simulations with the Heisenberg hamiltonian. The uniform 
ferromagnetic phase is obtained with the DFT exchange interactions at room temperature for the 
impurities concentration of 7 %. With the DFT-I-U exchange interactions, the ferromagnetic clusters 
form above room temperature in MC samples containing only 3 % Re. 
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I. INTRODUCTION 

Silicon doped with rhenium belongs to the widely stu- 
died class of materials called diluted magnetic semicon- 
ductors (DMS). In this class, the host materials belong to 
the group III-V, II- VI, and IV semiconductors or transi- 
tion metal oxides (DM0). Usually the dopants are the 
3rd-row d-elements, with Mn, Fe and Cr to be most 
frequently published. [1-6] Interestingly, the 5th-row el- 
ements, such as Re, show different properties, especially 
stronger hybridization with the host, and also deserve 
investigation. [7] DMS attract much attention, experi- 
mental and theoretical, because of their promissing ap- 
plication in spintronics.[8-10] To be used in fabrication 
of electronic devices, materials must preserve their prop- 
erties in conditions much above room temperature. To 
date, a very few reports on accessing room temperature 
has been published, and remain still controversial be- 
cause Curie temperature depends a lot on the sample 
preparation. [11-14] 

Amorphous silicon doped with rhenium is ferromag- 
netic at room temparature, but samples loose this prop- 
erty after a few months. [15] 

Basic theoretical information about single rhenium im- 
purities in crystallic silicon, in various crystal positions 
and charged states, have been recently reported. [7] How- 
ever, ferromagnetism in disordered materials is a com- 
plex phenomenon, where single-site impurity character- 
istics is as important as interactions between the mag- 
netic centers, and also their clustering. [4, 6, 16-19] There- 
fore, here, a one step closer to the experiment is done, 
and the work focus on magnetic interactions between two 
Rc impurities embeded in the silicon substitutional sites. 
The calculations are performed within the density func- 
tional theory. [20] Simultaneously, the effect of Coulomb 
interactions is investigated by means of the DFT-I-U 
scheme, [21] since it is well known that the strong interac- 
tions within the d-shell favour Hund's rules, change the 
local magnetizations and the binding between impurities, 
and they strongly influence also the anisotropy of inter- 



actions, as well as complex observables like the critical 
temperature. [6, 22] In the next step, the exchange inter- 
actions, obtained within the DFT and the DFT-I-U meth- 
ods, are used in Monte Carlo simulations of the Heisen- 
berg model and the Curie temperature is estimated. 

The main findings concern the geometric, magnetic 
and statistical effects. In physics of the short-distance 
interacting pairs, the local geometry plays a very impor- 
tant role, since two rhenium atoms attract each other and 
form the chemical bonds, which in turn affect very much 
magnetizations. Therefore, these cells are relaxed in the 
calculations. Further, the ferromagnetic (FM) and the 
antiferromagnetic (AF) energy difference is examined as 
a function of the Re-Re distance and orientation in the 
crystal. In both theoretical methods, the FM interactions 
are lower in the energy. Moreover, in the DFT-I-U scheme 
the magnetic couplings are stronger and very anisotropic. 
The mechanism of magnetic interactions between Re im- 
purities in silicon is the ferromagnetic superexchange. 
The formation energies calculated for all pairs lead to 
the close pair scenario with a very few single impurities. 
Clusters with more than two rhenium atoms are out of 
scopes of this theoretical studies, but in the future the 
agregates of the substitutional and the interstitial impu- 
rities will be investigated too. Finally, the transition from 
a magnetically disordered to an ordered phase has been 
studied by means of MC simulations. It was found that, 
this transition was different when the magnetic interac- 
tions calculated with the DFT or the DFT-f U methods 
were embeded to the model. Two main results are ob- 
tained: i) the uniform growth of the total magnetization 
with decreasing temperature using the pair interactions 
calculated within the DFT method, and ii) formation of 
magnetic clusters, when the DFT-I-U interactions were 
implanted into the MC simulations, similarly to the ef- 
fect of superparamagnetic blocking. [6] 

This paper is organized as follows: in section II the 
computational details are given, in section HI the results 
for close- and medium-distance pairs are collected, in sec- 
tion IV the long range effects and magnetic anisotropy 
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are described together with the formation energies of all 
Re-Re pairs, in section V the implication of calculated 
exchange interactions for the Curie temperature is dis- 
cussed, the summary is given in section VI. 

II. COMPUTATIONAL DETAILS 

The ab initio calculations of total energies were per- 
formed within the density functional theory, [20] us- 
ing the plane-wave pseudopotential code QUANTUM 
ESPRESSO. [23] The scalar relativistic pseudopotential for 
rhenium has been generated for the valence configuration 
6s2 5d5, whithin the ultrasof pseudopotential (USPP) 
scheme, [24] including the nonlinear core correction. [25] 
The DFT functional, within the generalized gradient ap- 
proximation (GGA) of Perdew-Wang (PW91) type,[26] 
has been employed for calculations of the electronic struc- 
ture of rhenium pairs in many large supercells. All cal- 
culations were repeated with the GGA+U method, [21] 
accounting for the strong correlations effects. The effec- 
tive U parameter (in fact U-J) has been estimated from 
the response method, [27] and its value of 2.3 eV is the 
same as in the previous work. [7] 

The plane-wave cutoffs for the energy of 30 Ry and 
for the density of 300 Ry are sufficient for the USPP 
without the semicore electrons in the valence calcula- 
tions. For the smearing at the Fermi siirfacc, the metal- 
lic broadening [2 8] of 0.01 Ry is used. The k-point mesh 
4x4x4 has been checked, to give the same results as 
the twice denser mesh. 

All supercells were constructed by repetition of the 
simple cubic 8-atom cell. The supercell sizes were: 
2x2x2 (64 atoms, the concentration of Re was x=3.125 
%), 3 X 2 X 2 (96 atoms, a;=2.083 %), 4 x 2 x 2 (128 
atoms, a;=1.5625 %), 3 x 3 x 2 (144 atoms, a;=1.3888 %) 
and 4x4x2 (256 atoms, x=0. 78125%). Above procedure 
has been applied in order to gather all possible directions 
with a minimal computational cost. The first Re atom 
was placed at the crystal site (0,0,0), and the second 
impurity was at the site (fc, /, m) in units of the silicon 
lattice constant divided by four. We did not account for 
interactions of the impurities implanted in the calculated 
cells with their periodic images (atoms from the neigh- 
boring cells), since for neutral cells these effects are not 
very pronounced. [29, 30] All calculations were performed 
at the lattice constant of 10.32 a.u., which was optimized 
with the applied pseudopotcntials. 

Monte Carlo (MC) simulations of the thermal average 
of magnetization (M) were performed using the Metropo- 
lis algorithm [31] with the Heisenberg hamiltonian: 

H = --'^JijSiSj. (1) 

ij 

The interaction constants Jij act between two impurities 
at sites i and j, and they are equal to the total energy 
difference in the antiferromagnetic and the ferromagnetic 



state (Eaf — Efm)- Numbers Si and Sj are magnetic 
moments (+1 or -1) at the doped sites. 

The cells for the simulations have been constructed 
from L X L X L replicas of the 8-atom simple cubic cell. 
The periodic boundary conditions were applied in MC su- 
percells. The critical temperatures (Tc) were estimated 
from the crossing of the 4th-order Binder cumulant[32] 
curves, defined as Ul{T) = 1 - {M^)/{3{MY), and cal- 
culated in cells of the size L =16, 20 and 24. The number 
of MC measurements was 5000 sweeps per impurity, and 
for the warming 500 sweeps per impurity. Some of MC 
simulations were repeated with higher number of sweeps 
per impurity, equal to 50000, for both the warming and 
the MC measurements, to check the convergence. The 
sites for the impurities have been chosen randomly in 
3D supercells according to an algorithm of the 3rd-order 
Sobolev's chains, [33] in order to obtain a quasi uniform 
random distribution. At the end, the average over 32 
random geometries was taken, to find the magnetization 
and the Ul curves. Simulations have been performed at 
the concentrations of 3, 5 and 7 % for the GGA pair in- 
teractions, and at lower concentrations of 2, 3 and 5 % 
for the GGA+U interactions (in this method the rhenium 
ions couple much stronger than in the GGA). 

III. RESULTS FOR CLOSE- AND 
MEDIUM-DISTANCE PAIRS 

In our previous work, [7] the lattice relaxation around 
the substitutional site was negligible, but for the impu- 
rities which interact on a short distance the geometry 
optimization is important. Therefore, the positions of 
atoms in the cells with pairs: 111, 220, 113, 331 and 333 
were optimized within the Newton- Raphson optimization 
scheme based on the Broyden-Fletcher-Goldfarb-Shanno 
(BFGS) algorithm. [34] It is found that the optimized Re- 
Re distances for the close pairs are shorter, and that 
for medium separations these distances almost do not 
change; the corresponding numbers are: 6.4 % for the 
111 pair, 23.6 % for 220, 3 % for 113, 0.2 % for 331 and 
0.05 % for 333. From now, the relaxed cells are denoted 
by "R", and results for the cells with considerable relax- 
ations, e.g. IIIR, 220R and 113R, will will be presented 
in detail. 

Table I collects the total and the absolute magnetiza- 
tions per Re and impurities separations in the cells with 
short- and medium-distance pairs. In the GGA method, 
the very close pairs, 111 and lllR, are nonmagnetic. The 
pair 220R also converges to the nonmagnetic state, if one 
starts from the AF state, while in the FM state it is mag- 
netic. The nonrelaxed pair 220 is magnetic in both the 
FM and the AF states, however the magnetization of the 
AF state is largely reduced. The same situation with 
even stronger reduction of local magnetic moment is in 
the relaxed cell 113R. In contrast, within the GGA+U 
method, almost all pairs are magnetic in both the FM 
and the AF state, except the pairs 111 and lllR which 
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TABLE 1. Total/absolute magnetizations (in /ib) in the cells 
in the FM state, and absolute magnetizations for the AF 
state, and the impurities separations Upair (in a.u.) 
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are nonmagnetic when one starts from the AF state. All 
magnetizations in the GGA+U method are higher than 
in the GGA method. This is usual behaviour, due to the 
larger exchange splitting whithin the e^ states of rhe- 
nium; see the density of states (DOS) of a single substi- 
tutional Re in the previous calculations. [7] 

In this point, it is worth to comment on nonmagnetic 
states in the close pairs: the substitutional Re in sili- 
con occupies halfly the minority spin Cg state, and this 
is different from the homovalent element Mn, where the 
crystal field splitting is weaker than the exchange split- 
ting and the \,2g states are at the Fermi level. [7] In the Mn 
doped silicon, the close pairs are more strongly magne- 
tized than the separate pairs, due to partial occupation of 
majority spin bands. When the minority spin bands are 
partially occupied, then the close interactions increase 
filling of these bands to more "closed-shell" scenario, and 
the spin polarization decreases. 

The density of states for the nonmagnetic solution of 
the ferromagnetic start for the 111 pair and for the FM 
phase of the 220 pair obtained within the GGA and 
GGA+U methods are presented in Figurel. The strong 
hybridization of 5d-Re states with the neighboring atoms 
is obvious for the 111 pair. In the case of 220 pair, the 
degeneracy of the e^, spin-down states is partially lifted 
due to the broken symmetry around the Re atom when 
the second impurity is placed very closely. 

IV. LONG-DISTANCE PAIRS AND MAGNETIC 
ANISOTROPY AND FORMATION ENERGIES 

This section is focused on general trends in magnetic 

interactions for all pairs, and on the formation energies 
with respect to the most bound pair lllR; Table II gath- 
ers all these numbers. 

The formation energy of all pairs (except 111, which 
however is only given for the comparison, to illustrate the 
relaxation effect) cost about 1 eV more than formation 
of the close lllR dimer. The fact, that the GGA+U 
formation energies are lower than those from the GGA, 



FIG. 1. The density of states (DOS) for the 111 and 220 pairs 
in the FM/nonmagnetic states obtained within the GGA and 
the GGA+U methods. The \,2g and Cg states for one of the 
impurities in the pair arc drawn on the gray color background 
of the 5d-states of both Re atoms. The Fermi level is marked 
by the vertical line. 

does not change the scenario of very close pairs. At this 
point, one could say, that it is not worth to deal with the 
rest of pairs because statistically there will be a very few 
of them in the sample. It would be reasonable to focus 
on clusters with more than two impurities and it will be 
a subject of the future work. Now, the conclusions from 
medium- and long-distance pairs will be drawn, since in 
the large sample they also exist and mediate magnetic 
interactions for a long distance, as will be discussed in 
connection with the critical temperature. 

The magnetic couplings Jint=EAF — Efm, obtained 
with the GGA method, decay quickly with the distance, 
and only the pairs in the 64 atom cell are strongly cou- 
pled, with the maximal exchange interaction found for 
the 220 pair. All pairs are coupled ferromagnetically 
(positive numbers Jmt) iu both applied here methods. 

As for the anisotropy of magnetic interactions, the 
GGA numbers do not show any visible dependence on 
the crystal direction. While among the GGA+U values, 
there is well pronounced enhancement of the coupling 
along the 110 crystal axis. This direction is prefered 
for the couplings in zinc-blende hosts, due to the zygzak 
chains, which mediate magnetic interactions via a small 
polarization of the host atoms neighboring to impurities. 
These polarizations are visible in the Lowdin population 
analysis, [35] and amount to about 0.01+0.02 of electron 
charge per the Si site. 

The range of magnetic interactions in the GGA+U 
method is long. There are quite large J,;„/ numbers 
along 110-direction for the distance over 20 a.u.; these 
are for instance numbers 23 meV for the 660 pair, and 
about 6 meV for the 880 pair in the 256-atoms cell. 
At long-distances, even small couplings play a very cru- 
cial role in estimation of the Curie temperature, because 
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TABLE II. Interactions Jint (in meV) between rhenium pairs. 
We assumed Jint being negligible for pairs: 755, 775, 777, 862, 

866, 884, 888 (256 atoms cell). Ef is the formation energy (in 
meV) of the lower-energy magnetic phase for the given pair, 
with respect to the formation energy of the dimer 111_R. The 
distance Rpair (in a.u.) between Re-Re, and the coordination 
number Npair of given distance impurities within the shell 
around the impurity ion at (0,0,0) are also given. 

they permit a percolation between impurities at small 
concentrations. [36] For instance, it is known, that to get 

FM-ordcr at concentration of 7 %, the exchange interac- 
tions must extend to at least five coordination shells. [6] 
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TABLE III. Curie temperatures (in K) at various concentra- 
tions X obtained from the MC simulations by: 1) the magne- 
tization curves crossing ^^'^('^"i and 2) the cumulant curves 
crossing T^'^^'^^\ of averaged samples, using the GGA inter- 
actions Jint from Table II. Tm>o is the temperature range 
at which the magnetization in studied samples starts to rise 
from zero. 

V. CRITICAL TEMPERATURES 

In this section, the discussion of the results of Monte 
Carlo simulations with the Heisenberg model is pre- 
sented. The simulations were performed on a number 
of supercells with different sizes, and with various im- 
purities distributions, and for a few concentrations. In 
the hamiltonian, the Re- Re interactions Jint, calculated 
within the GGA and the GGA+U schemes (see Table II) 
have been used. 

For the interactions Jint obtained within the GGA 
method, which are rather not anisotropic in space, all 
samples start to magnetize at some temperature, and 
this magnetization saturates to high values of 0.7-^0.95 
per impurity, depending on a sample. The starting tem- 
perature for a visible rise of magnetization depends quite 
strongly on a geometry of impurities distribution. The 
Curie temperatures for the ferromagnetic order were esti- 
mated in two ways: (1) from the crossing of the magneti- 
zation curves, and (2) from the crossing of the 4th-order 
cumulant curves. [32] These curves have been calculated 
for different supercell sizes, and each of these curves has 
been averaged over geometries of impurities distribution 
- as it is described in section with the computational 
details. 

The Curie temperatures, calculated with the aforemen- 
tioned procedure, are presented in Table III. The temper- 
atures obtained from the magnetization-curves crossing 
are a bit lower than Tc obtained from the cumulant- 
curves crossing. The last column shows the temperature 
at which the magnetization in samples starts to rise from 
zero {T^Myo)' this temperature is not averaged over 
the impurities distributions, but it is given as a range to 
which fall the results from all used geometric samples. 
This quantity shows how much the theoretical results for 
given concentration differ just due to different impurities 
distribution. It gives also an argument in the explana- 
tion why the experimental Tc, published for the same 
material and the concentration, vary a lot depending on 
samples preparation. 

In contrast to simulations with the GGA interactions, 
the GGA-I-U interactions implanted to the Heisenberg 
hamiltonian usually do not lead to the total and uniform 
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TABLE IV. Approximate temperature ranges T^fm (in K) 
at which the ferromagnetic cluster order starts to form in 

the MC supercells containing impurities at concentration x, 
and using the interactions Jint calculated by means of the 
GGA+U method (from Table II). 

magnetization in the supcrccll. Independently on geo- 
metric parameters used in these simulations, it is found 
negligible or very low magnetization of samples (0.1-=-0.3 
jiB per impurity) at most simulation temperatures. In- 
terestingly, instead of homogeneous magnetization within 
a supcrccll, the ferromagnetic clusters form. In some 
samples, there was a transition temperature at which it 
was found a large total magnetization and below this 
temperature the magnetic clusters formed. In other sam- 
ples, there is a smooth transition between magnetically 
disordered phase and the phase with cluster structure. 
The magnetization of local clusters is also not full, but 
grows with the decreasing temperature. One can only 
approximate the transition temperature by watching the 
pictures for different samples. 

Figure 2 presents magnetization of one chosen geomet- 
ric sample, calculated using the interactions from the 
GGA-I-U method, with impurities concentration of 5 %, 
at five temperatures: 116 K, 464 K, 696 K, 812 K and 
928 K. The sample size is LxLxLxS with L=20, and it 
is divided into a few boxes, and each of these boxes is 
of the size 4x4x4x8; thus, the whole supecell is divided 
into 5x5x5 pieces. Whithin the each small box, the mag- 
netization over impurities is averaged, and it is done over 
the spins obtained in the last MC measurement sweep. 
The 3D supcrccll is presented in the slices, which cut 
parallel to the XY-plane, and these cut slices arc placed 
in Figure 2 in columns. The averaged magnetization per 
impurity can take the value from -1 to 1, therefore the 
absolute value of box-magnetization is marked by a circle 
with radius proportional to magnetization, and spin di- 
rection (up and down) is depicted by the white and black 
color. 

One can see in Figure 2 that the FM-ordcr transition, 
in presented sample, occurs between 820 and 900 K. In 
some other samples (not presented in Figure 2), at the 
same concctration of 5 %, the estimation of T^^ is not 
so easy because one cannot find a temperature at which 
the total sample magnetize; only the magnetic clusters 
slowly grow. 

It is important to mention, that this scenario with 
magnetic clusters does not change with a very large num- 
ber of measurement sweeps, of 50 000 per impurity, and 
with the same number of warming sweeps in the MC sim- 
ulations. 



In Table IV, there are listed the rouly estimated tem- 
peratures at which the cluster FM-ordcr starts to form, 
'^¥mi f'^'" concentrations: 2, 3, and 5 %. The range of 
these temperatures is wide, because the magnetic tran- 
sition depends on impurities distribution in the sample, 
like it was the case with the GGA critical temperatures. 
At the concentrations above 4-5 %, all samples form the 
magnetic cluster structure above room temperature. 

It is known about the percolation that, if impu- 
rities distribution is of cluster type and clusters are 
large, then a superparamagnetic blocking phenomenon 
becomes important. [6, 37] In the simulations performed 
in this work, the distributions are uniform nctherveless 
the superparamagnetic blocking occurs, and the reason 
for this is the magnetic coupling of pairs which acts at 
a very long range. When during the MC simulations 
the magnetic clusters form, their internal coupling is so 
strong, that it costs much energy to flip the whole cluster. 
Finally, the magnetic relaxation time is very long. Thus, 
the total magnetization, in some of the samples with the 
cluster structure, is close to zero over all temperatures, 
even after 50 000 MC measurements per impurity. 

VI. SUMMARY 

In this work, the rhenium pairs cmbcded in silicon in 
substitutional sites were calculated by means of the GGA 
and the GGA-I-U methods. The formation energy of the 
closest pair is roughly 1 eV lower than those energies 
for the rest of double combinations. This supports the 
cluster scenario with a very few single impurities and 
medium-distance pairs. It is a common behaviour of im- 
purities in various DMS.[16 19] 

The calculated magnetic couplings iint=^AF-£'FM are 
positive for all pairs, and very long range. The GGA 
couplings arc not anisotropic with respect to the crys- 
tal axes, in contrast the GGA+U interactions are much 
stronger along the 110 direction. At small distances, the 
Re-Re pairs become nonmagnetic for both the AF and 
FM phases, or only for the AF phase. At higher concen- 
trations this effect is even stronger, but less pronounced 
in the GGA+U method than in the GGA scheme. 

The mechanism which drives this DMS magnetic is the 
ferromagnetic superexchange, with high critical temper- 
atures excessing room temperature, even for low concen- 
trations of impurities of 7 % for the GGA and 3 % for the 
GGA+U. Magnetizations of Monte Carlo samples below 
critical temperatures arc uniform for the GGA, and dif- 
ferent with magnetic clusters structure for the GGA+U 
scheme. Quite large and very extended magnetic cou- 
plings support percolation in this DMS, and give hope to 
obtain still high Curie temperature when the impurities 
clusters are included in simulations. 

There is also plaussible the existence of superparam- 
agnetic blocking phenomenon below the Curie temper- 
ature. But this scenario must be confirmed by further 
studies with clusters of the substitutionals and the inter- 
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stitials together, and also employing a method for better 
treatment of the self-interaction effect in the host (sil- 
icon), namely the DFT with the pseudopotential self- 
interaction correction (pSIC).[38] 

VII. ACKNOWLEDGMENTS 

I thank Prof. R. R. Galazka for showing me the exper- 
imental results and to Andrzej Fleszar for many philo- 



sophical discussions. The calculations have been per- 
formed at the Leibniz Supercomputing Centre in Mu- 
nich. The work was supported by the European Founds 
for Regional Development within the SICMAT Project 
(Contract No. UDA-POIG.01.03.01-14-155/09). 



[1] J. K. Furdyna, J. Appl. Phys. 64(4), R29 (1988). 

[2] J. K. Furdyna, Semimagnetic Semiconductors and Di- 
luted Magnetic Semiconductors, ed M. Averous and 
M. Balkanski (New York: Plenum) 1991. 

[3] H. Ohno, Science 281, 951 (1998). 

[4] T. Dietl, Semicond. Sci. Technol. 17, 377 (2002). 

[5] S. Sanvito, G. Theurich and N. A. Hill, Journal of Super- 
conductivity and Novel Magnetic Materials 15, 85 (2002) 

[6] K. Sato, L. Bergqvist, J. Kudrnovsky, P. H. Ded- 
erichs, O. Eriksson, 1. Turok, B. Sanyal, G. Bouzcrar, 
H. Katayama-Yoshida, V. A. Dinh, T. Fukushima, 
H. Kizaki, R. Zcllcr, Rev. Mod. Phys. 82, 1633 (2010). 

[7] M. Wierzbowska and A. Fleszar, Phys. Rev. B 83, 184418 
(2011). 

[8] R. R. Galazka and J. Kossut, Lecture Notes in Physics 

(Springer, Berlin, 1982), Vol. 132, p. 245. 
[9] H. Ohno, Nature Materials 9, 952 (2010). 

[10] T. Dietl, Nature Materials 9, 965 (2010). 

[11] M.L. Reed, N.A. El-Masry H.H. Stadelmaier, M.K. Rit- 
ums, M.J. Reed, C.A. Parker, J.C. Roberts, and S.M. Be- 
dair, Appl. Phys. Lett. 79, 3473 (2001). 

[12] M. Bolduc, C. Awo-Affouda, A. StoUenwerk, 
M. B. Huang, F. G. Ramos, G. Agnello, and V. P. La- 
Bella, Phys. Rev. B 71, 033302 (2005). 

[13] Yuji Matsumoto, Mokoto Murakami, Tomoji Shono, Tet- 
suja Hasegawa, Tomoteru fukumura, Masashi Kawasaki, 
Parhat Ahmet, Toyohiro Chikyow, Sliin-ya Koshihara, 
and Hideomi Koinuma, Science 291, 854 (2001). 

[14] J. M. D. Coey, M. Viret and S. vou Mohiar, Adv. Phys. 
48, 167 (1999); Y. Tokura, Colossal Magnetoresistive Ox- 
ides, ed. Y. Tokura (London: Gordon and Breach). 

[15] R. R. Galazka et al., to be published. 

[16] T. Hynninen, H. Raebiger, J. von Boehm, Phys. Rev. B 
75, 125208 (2007). 

[17] N. Gonzalez Szwacki, J. A. Majewski, and T. Dietl, 
Phys. Rev. B, 83, 184417 (2011). 

[18] A. Navarro-Quczada, N. Gonzalez Szwacki, W. Stefanow- 
icz, Tian Li, A. Grois, T. Devillers, R. Jakiela, B. Faina, 
J. A. Majewski, M. Sawicki, T. Dietl, and A. Bonanni, 
arXiv: 1107.4901 

[19] V. N. Monshov, V. V. Tugushev, S. Caprara, and 
E. V. Chulkov, Phys. Rev. B, 83, 035201 (2011). 

[20] P. Hohenborg and W. Kohn, Phys. Rev. 136, B864 

(1964) ; W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 

(1965) . 

[21] V.I. Anisimov, J. Zaanen, and O.K. Andersen, 
Phys. Rev. B 44, 943 (1991); V.I. Anisimov, F. Aryaseti- 
awan, and A.I. Lichtenstein, J. Phys.: Condens. Matter 



9, 767 (1997). 

[22] K. Sato, P. H. Dederichs, and H. KatayamarYoshida, 
Physica B 376-377, 639 (2006); M. Toyoda, H. Akai, 
K. Sato, and H. Katayama-Yoshida, Physica B 376-377, 
647 (2006). 

[23] P. Giannozzi ct al., J. Phys. Condens. Matter, 21, 395502 
(2009). 

[24] D. Vanderbilt, Phys. Rev. B 41, R7892 (1990). 

[25] S. G. Louie, S. Proyen, and M. L. Cohen, Phys. Rev. B 

26, 1738 (1982). 
[26] J. P. Perdcw and Y. Wang, Phys. Rev. B 45, 13244 

(1992). 

[27] M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 
035105 (2005). 

[28] N. D. Mermin, Phys. Rev. A 137, 1441 (1965); 

M. J. Gillan, J. Phys. Condens. Matter 1, 689 (1989). 
[29] G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 

(1995) . 

[30] M. I. J. Probert and M. C. Payne, Phys. Rev. B 67, 
75204 (2003). 

[31] N. Metropolis, A. W. Rosenbluth, M. n. Rosenbluth, 
A. H. Teller and E. Teller, J. Chem. Phys. 21, 1087 

(1953). 

[32] K. Binder, Phys. Rev. Lett. 47, 693 (1981). 

[33] Numerical Recipies in Fortran, edited by W. H. Press, 
S. A. Teukolsky, W. T. Vetterling, B. P. Flannery - sec. 
cd. (Cambridge University Press 1986 1992), p. 299. 

[34] Roger Fletcher, Practical Methods of Optimization, John 
Wiley and Sons, Chichester, 2nd edition, 1987. 

[35] A. Szabo and N. S. Ostlund, Modern Quantum Chem- 
istry. Introduction to Advanced Electronic Structure The- 
ory., Dover Publications INC., Mineola, New York 

(1996) . 

[36] A. Kaminski and S. Das Sarma, Phys. Rev. Lett. 88, 
247202 (2002). 

[37] A. Aharoni, Introduction to the Theory of Ferromag- 

netism, Oxford University Press, Oxford. 
[38] A. Filippetti and V. Fiorentini, Eur. Phys. J. B 71, 139 

(2009). 



T=116K 




T=464K 



•• 






r ^ 








r ^ 




• 






• 


• 


• 






• 


• 


• 







T=696K 





■ 


• 


• 




■ 


• 


• 




□ 


• 


• 




■ 


• 


• 




■ 


• 


• 



T=812K 



T=928K 

^•1 ■ 





• 




• 
























• 




□ 




• 




■ 








■ 




W 





□ 


• 


• 




D 


• 


• 




□ 


• 


• 


3: 


D 


• 


• 


ri 


• 


■ 

















• 


nr 






• 








• 











nn 


■ 


• 


• 




H 


• 
• 


• 

• 




B 


• 


• 




■ 


• 


• 



□ 




■ 




• 


D 




■ 


• 


• 


□ 




D 


• 


• 


□ 




□ 


• 








□ 


m 


■ 






B 


• 




nr ^ 


■ 


• 


• 




■ 


• 


• 




■ 


• 


• 


r-ir ^ 


□ 


• 


• 



FIG. 2. Ferromagnetic clusters in one chosen MC sample, 
which contains 5 % Re. Simulations were performed with the 
GGA+U interactions from Table II and at five temperatures: 
116 K, 464 K, 696 K, 812 K, 928 K. Each column presents five 
cuts through the sample, which are perpendicular to the xy- 



